Spatial Satistics
Fernando García
1/10/2020
Introducción
Referencia de la imagen: “The TibA Adhesin/Invasin from Enterotoxigenic Escherichia coli Is Self Recognizing and Induces Bacterial Aggregation and Biofilm Formation”
Asumimos que los puntos pueden tener cualquier localización en la ventana, en nuestro caso no hay ninguna limitación al respecto. El objetivo es estimar los parámetros de distribución de las diferentes bacterias.
En los datos tenemos dos tipos de poblaciones bacterianas (marks = 2).
#install.packages("imager")
#install.packages('spatstat')
library(imager)
library(spatstat)
library(cowplot)
library(ggplot2)image <- load.image('./Bacterial_interaction.png')
plot(image, main = "Original image")
image <- resize(image,round(width(image)/4),round(height(image)/4))
plot(image,main="Low resolution") 
image## Image. Width: 112 pix Height: 112 pix Depth: 1 Colour channels: 4Pasar la imagen a un dataFrame.
image.df <- as.data.frame(image)
dim(image.df)## [1] 50176 4Selección de los puntos con mayor intensidad de acuerdo al color verde y rojo (RGB). Con esta técnica algunos puntos de la imagen original pueden quedar excluidos. Se han seleccionado aquellos que en el correspondiente canal tengan más de un 50% de intensidad y el resto de canales menos del 0.1. Se ha realizado un filtrado de los valores que en el canal 3 (azul) tiene una intensidad superior a 0.3
image.df$quality <- image.df$cc
image.df[image.df$cc == 3 & image.df$value > 0.2, ]$quality = 0spatial.points <- data.frame(x = image.df[image.df$value > 0.5 & (image.df$cc == 2 | image.df$cc == 1) & image.df$quality != 0, ]$x,
y = image.df[image.df$value > 0.5 & (image.df$cc == 2 | image.df$cc == 1) & image.df$quality != 0, ]$y,
strain = image.df[image.df$value > 0.5 & (image.df$cc == 2 | image.df$cc == 1) & image.df$quality != 0, ]$cc)p <- ggplot(spatial.points, aes(x=x, y=y)) +
geom_point(aes(color = as.factor(strain))) + theme_minimal() + theme(legend.title = element_blank()) +
scale_color_manual(values = c("#FF0000", "#04FF00"))
p
Carga de los puntos con una ventana de dimensiones 112 x 112
pattern <- ppp(spatial.points$x, spatial.points$y, c(0, 112), c(0, 112))
plot(pattern, main = NULL)
summary(pattern)## Planar point pattern: 716 points
## Average intensity 0.05707908 points per square unit
##
## *Pattern contains duplicated points*
##
## Coordinates are integers
## i.e. rounded to the nearest unit
##
## Window: rectangle = [0, 112] x [0, 112] units
## Window area = 12544 square unitspattern <- rescale(pattern, 10) # Rescale to lower distance
summary(pattern)## Planar point pattern: 716 points
## Average intensity 5.707908 points per square unit
##
## *Pattern contains duplicated points*
##
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 units
##
## Window: rectangle = [0, 11.2] x [0, 11.2] units
## Window area = 125.44 square unitsIntensidad
intensity(pattern)## [1] 5.707908Error estandar
sqrt(intensity(pattern)/area(Window(pattern)))## [1] 0.2133145Quadrat count
Q <- quadratcount(pattern, nx=6, ny=6)
Q## x
## y [0,1.87) [1.87,3.73) [3.73,5.6) [5.6,7.47) [7.47,9.33)
## [9.33,11.2] 17 32 24 21 6
## [7.47,9.33) 21 29 15 18 12
## [5.6,7.47) 9 12 27 9 10
## [3.73,5.6) 28 13 19 21 31
## [1.87,3.73) 18 17 30 10 10
## [0,1.87) 21 27 21 21 25
## x
## y [9.33,11.2]
## [9.33,11.2] 33
## [7.47,9.33) 21
## [5.6,7.47) 17
## [3.73,5.6) 18
## [1.87,3.73) 24
## [0,1.87) 29plot(intensity(Q, image = TRUE), main = "Intensity per quadrant GFP + RFP")
La hipótesis nula es que el patrón se distribuye de forma homogénea. En nuestro caso podemos rechazar dicha hipótesis de forma que aceptamos que nuestros datos no son homogéneos.
quadrat.test(pattern, nx=6, ny=6)##
## Chi-squared test of CSR using quadrat counts
##
## data: pattern
## X2 = 95.609, df = 35, p-value = 3.13e-07
## alternative hypothesis: two.sided
##
## Quadrats: 6 by 6 grid of tilesVolviendo a realizar el test con un número menor de cuadrantes tenemos que se acepta la hipótesis nula.
quadrat.test(pattern, nx = 5, ny = 5)##
## Chi-squared test of CSR using quadrat counts
##
## data: pattern
## X2 = 47.617, df = 24, p-value = 0.005628
## alternative hypothesis: two.sided
##
## Quadrats: 5 by 5 grid of tilesquadrat.test(pattern, nx = 4, ny = 4)##
## Chi-squared test of CSR using quadrat counts
##
## data: pattern
## X2 = 21.654, df = 15, p-value = 0.2344
## alternative hypothesis: two.sided
##
## Quadrats: 4 by 4 grid of tilesKernel density
Se trata de una estimación no paramétrica de la función de intensidad.
d <- density(pattern, edge = TRUE, kernel="gaussian")
plot(d, main = "Gaussian kernel density stimation GFP + RFP")
Dependencia entre puntos
Este modelo asume que el proceso es estacionario (proceso estocástico con una distribución de probabilidad relativamente constante).
rescaled.pattern <- rescale(pattern, 10)
fryplot(rescaled.pattern, main = "Fry GFP + RFP")
La función empírica K asume asume de nuevo un proceso estacionario.
plot(Kest(pattern), main="K function")
Cada función corresponde a diferentes correcciones de los ejes: * theo Kpois(r) theoretical Poisson K(r)
* border Kbord(r) border-corrected estimate of K(r)
* trans Ktrans(r) translation-corrected estimate of K(r) * iso Kiso(r) isotropic-corrected estimate of K(r)
Transofrmando la función K a un alinea recta mediante la función L de Besag obtenemos la siguiente gráfica que nos permite apreciar mejor las desviaciones. La función aplica las mismas correcciones que aplicaba antes.
plot(Lest(pattern), main = "Besag’s L-function")
La función de correlacción por pares (pair correlation function) considera en el numerador la probabilidad de observar un par de puntos separados por una distancia r, dividida por esa misma probabilidad para un proceso estocástico. En este caso se puede ver que hay un claro patrón de clustering dentro de una pequeña distancia. Esto puede ser indicativo de posibles agrupaciones entre un número reducido de bacterias.
plot(pcf(pattern), main = NULL)
Se han seleccionado 99 simulaciones, de esta forma, la hipótesis nula es rechazada con un valor de alfa = 0.01.
a <- capture.output(plot(envelope(pattern, Kest, nsim = 99), main = NULL, legend = FALSE))
a <- capture.output(plot(envelope(pattern, Lest, nsim = 99), main = NULL, legend = FALSE))
Espaciado entre los puntos
La información referente al espaciado entre los puntos proporciona información complementaria al estudio de la correlación realizado anteriormente. La función empty-space (función F) asumiendo de nuevo un estado estacionario corresponde a la función de distribución empírica de las distancias observaas respecto al espacio vacio de m localizaciones. En este caso tenemos una distancia acumulativa. De nuevo la librería spatstats aplica diferentes correcciones.
plot(Fest(pattern), main = "F-function (point-to-event)")
a <- capture.output(plot(envelope(pattern, Fest, nsim = 99), main = NULL, legend = FALSE))
La función G (Nearest neighbour distances o event-to-event) de nuevo asumiendo un estado estacionario computa la probabildiad de que un punto tenga una distancia menor, mayor o igual a un proceso estacionario.
plot(Gest(pattern), main = "G-function (event-to-event)")
a <- capture.output(plot(envelope(pattern, Gest, nsim = 99), main = NULL, legend = FALSE))
La función J es una combinación de las funciones F y G.
plot(Jest(pattern), main = "J-function")
Considerando el patrón no homogéneo
En este caso suponiendo que los datos no tengan un patrón homogéneo de intensidad.
plot(Linhom(pattern), main = "L-function")
plot(Jinhom(pattern), main = "J-function")
plot(Ginhom(pattern), main = "G-function")
plot(Finhom(pattern), main = "F-function")
a <- capture.output(plot(envelope(pattern, Linhom, nsim = 99), main = "L-function", legend = FALSE))
a <- capture.output(plot(envelope(pattern, Jinhom, nsim = 99), main = "J-function", legend = FALSE))
a <- capture.output(plot(envelope(pattern, Ginhom, nsim = 99), main = "G-function", legend = FALSE))
a <- capture.output(plot(envelope(pattern, Finhom, nsim = 99), main = "F-function", legend = FALSE))
Análisis del patrón de las diferentes cepas bacterianas independientemenete
La metodología seguida ha sido la misma pero aplicada de forma independiente a cada una de las bacterias.
green.data <- data.frame(x = image.df[image.df$value > 0.5 & image.df$cc == 2 & image.df$quality != 0, ]$x,
y = image.df[image.df$value > 0.5 & image.df$cc == 2 & image.df$quality != 0, ]$y,
strain = image.df[image.df$value > 0.5 & image.df$cc == 2 & image.df$quality != 0, ]$cc)
red.data <- data.frame(x = image.df[image.df$value > 0.5 & image.df$cc == 1 & image.df$quality != 0, ]$x,
y = image.df[image.df$value > 0.5 & image.df$cc == 1 & image.df$quality != 0, ]$y,
strain = image.df[image.df$value > 0.5 & image.df$cc == 1 & image.df$quality != 0, ]$cc)p.green <- ggplot(green.data, aes(x=x, y=y)) +
geom_point(aes(color = as.factor(strain))) + theme_minimal() + theme(legend.title = element_blank()) +
scale_color_manual(values = c("#04FF00"))
p.red <- ggplot(red.data, aes(x=x, y=y)) +
geom_point(aes(color = as.factor(strain))) + theme_minimal() + theme(legend.title = element_blank()) +
scale_color_manual(values = c("#FF0000"))
plot_grid(p, p.green, p.red, nrow = 1, labels = c("A", "B", "C"))
Análisis de la intensidad
pattern.green <- ppp(green.data$x, green.data$y, c(0, 112), c(0, 112))
pattern.green <- rescale(pattern.green, 10)
pattern.red <- ppp(red.data$x, red.data$y, c(0, 112), c(0, 112))
pattern.red <- rescale(pattern.red, 10) intensity(pattern.green)## [1] 5.030293sqrt(intensity(pattern.green)/area(Window(pattern.green)))## [1] 0.2002528intensity(pattern.red)## [1] 0.6776148sqrt(intensity(pattern.red)/area(Window(pattern.red)))## [1] 0.07349764Q.green <- quadratcount(pattern.green, nx=6, ny=6)
Q.red <- quadratcount(pattern.red, nx=6, ny=6)
quadrat.test(pattern.green, nx = 6, ny = 6)##
## Chi-squared test of CSR using quadrat counts
##
## data: pattern.green
## X2 = 99.441, df = 35, p-value = 8.652e-08
## alternative hypothesis: two.sided
##
## Quadrats: 6 by 6 grid of tilesquadrat.test(pattern.red, nx = 4, ny = 4)##
## Chi-squared test of CSR using quadrat counts
##
## data: pattern.red
## X2 = 38.671, df = 15, p-value = 0.001435
## alternative hypothesis: two.sided
##
## Quadrats: 4 by 4 grid of tilesplot(intensity(Q.green, image = TRUE), main = "Intensity per quadrant in GFP")
plot(intensity(Q.red, image = TRUE), main = "Intensity per quadrant in RFP")
plot(density(pattern.green, edge = TRUE, kernel="gaussian"),
main = "Gaussian kernel density stimation in GFP")
plot(density(pattern.red, edge = TRUE, kernel="gaussian"),
main = "Gaussian kernel density stimation in RFP")
fryplot(pattern.green, main = "Fry GFP")
fryplot(pattern.red, main = "Fry RFP")
Análisis de la dependencia
plot(Kest(pattern.green), main="K function GFP")
plot(Lest(pattern.green), main="L function GFP")
plot(Kest(pattern.red), main="K function RFP")
plot(Lest(pattern.red), main="L function RFP")
plot(pcf(pattern.green), main = "GFP")
plot(pcf(pattern.red), main = "RFP")
a <- capture.output(plot(envelope(pattern.green, Kest, nsim = 99), main = "K-function GFP", legend = FALSE))
a <- capture.output(plot(envelope(pattern.green, Lest, nsim = 99), main = "L-function GFP", legend = FALSE))
a <- capture.output(plot(envelope(pattern.red, Kest, nsim = 99), main = "K-function RFP", legend = FALSE))
a <- capture.output(plot(envelope(pattern.red, Lest, nsim = 99), main = "L-function RFP", legend = FALSE))
Análisis de la distribución espacial
plot(Fest(pattern.green), main = "F-function GFP")
plot(Gest(pattern.green), main = "G-function GFP")
plot(Jest(pattern.green), main = "J-function GFP")
plot(Fest(pattern.red), main = "F-function RFP")
plot(Gest(pattern.red), main = "G-function RFP")
plot(Jest(pattern.red), main = "J-function RFP")
a <- capture.output(plot(envelope(pattern.green, Fest, nsim = 99), main = "F-function GFP"))
a <- capture.output(plot(envelope(pattern.green, Gest, nsim = 99), main = "G-function GFP"))
a <- capture.output(plot(envelope(pattern.green, Jest, nsim = 99), main = "J-function GFP"))
a <- capture.output(plot(envelope(pattern.red, Fest, nsim = 99), main = "F-function RFP"))
a <- capture.output(plot(envelope(pattern.red, Gest, nsim = 99), main = "G-function RFP"))
a <- capture.output(plot(envelope(pattern.red, Jest, nsim = 99), main = "J-function RFP"))